The data consists of vegetation % cover by functional group from across CONUS (from AIM, FIA, LANDFIRE, and RAP), as well as climate variables from DayMet, which have been aggregated into mean interannual conditions accross multiple temporal windows.

Dependencies

User defined parameters

print(params)
## $run
## [1] TRUE
## 
## $test_run
## [1] FALSE
## 
## $save_figs
## [1] FALSE
## 
## $ecoregion
## [1] "shrubGrass"
## 
## $response
## [1] "TotalHerbaceousCover"
# set to true if want to run for a limited number of rows (i.e. for code testing)
test_run <- params$test_run
save_figs <- params$save_figs
response <- params$response
fit_sample <- TRUE # fit model to a sample of the data
n_train <- 5e4 # sample size of the training data
n_test <- 1e6 # sample size of the testing data (if this is too big the decile dotplot code throws memory errors)


run <- params$run
# set option so resampled dataset created here reproduces earlier runs of this code with dplyr 1.0.10
source("../../../Functions/glmTransformsIterates.R")
source("../../../Functions/transformPreds.R")
source("../../../Functions/StepBeta_mine.R")
#source("src/fig_params.R")
#source("src/modeling_functions.R")
 
library(ggspatial)
library(terra)
library(tidyterra)
library(sf)
library(caret)
library(tidyverse)
library(GGally) # for ggpairs()
library(pdp) # for partial dependence plots
library(gridExtra)
library(knitr)
library(patchwork) # for figure insets etc. 
library(ggtext)
library(StepBeta)
theme_set(theme_classic())
library(here)
library(rsample)
library(kableExtra)
library(glmnet)
library(USA.state.boundaries)

read in data

Data compiled in the prepDataForModels.R script

here::i_am("Analysis/VegComposition/ModelFitting/02_ModelFitting.Rmd")
modDat <- readRDS( here("Data_processed", "CoverData", "DataForModels_spatiallyAveraged_withSoils_noSf.rds"))
## there are some values of the annual wet degree days 5th percentile that have -Inf?? change to lowest value for now? 
modDat[is.infinite(modDat$annWetDegDays_5percentile_3yrAnom), "annWetDegDays_5percentile_3yrAnom"] <- -47.8
## same, but for annual water deficit 95th percentile 
modDat[is.infinite(modDat$annWaterDeficit_95percentile_3yrAnom), "annWaterDeficit_95percentile_3yrAnom"] <- -600

# # Convert total cover variables into proportions (for later use in beta regression models) ... proportions are already scaled from zero to 1
# modDat <- modDat %>%
#   mutate(TotalTreeCover = TotalTreeCover/100,
#          CAMCover = CAMCover/100,
#          TotalHerbaceousCover = TotalHerbaceousCover/100,
#          BareGroundCover = BareGroundCover/100,
#          ShrubCover = ShrubCover/100
#          )
# For all response variables, make sure there are no 0s add or subtract .0001 from each, since the Gamma model framework can't handle that
modDat[modDat$TotalTreeCover == 0 & !is.na(modDat$TotalTreeCover), "TotalTreeCover"] <- 0.0001
modDat[modDat$CAMCover == 0 & !is.na(modDat$CAMCover), "CAMCover"] <- 0.0001
modDat[modDat$TotalHerbaceousCover == 0  & !is.na(modDat$TotalHerbaceousCover), "TotalHerbaceousCover"] <- 0.0001
modDat[modDat$BareGroundCover == 0 & !is.na(modDat$BareGroundCover), "BareGroundCover"] <- 0.0001
modDat[modDat$ShrubCover == 0 & !is.na(modDat$ShrubCover), "ShrubCover"] <- 0.0001
modDat[modDat$BroadleavedTreeCover_prop == 0 & !is.na(modDat$BroadleavedTreeCover_prop), "BroadleavedTreeCover_prop"] <- 0.0001
modDat[modDat$NeedleLeavedTreeCover_prop == 0 & !is.na(modDat$NeedleLeavedTreeCover_prop), "NeedleLeavedTreeCover_prop"] <- 0.0001
modDat[modDat$C4Cover_prop == 0 & !is.na(modDat$C4Cover_prop), "C4Cover_prop"] <- 0.0001
modDat[modDat$C3Cover_prop == 0 & !is.na(modDat$C3Cover_prop), "C3Cover_prop"] <- 0.0001
modDat[modDat$ForbCover_prop == 0 & !is.na(modDat$ForbCover_prop), "ForbCover_prop"] <- 0.0001
# 
# modDat[modDat$TotalTreeCover ==1& !is.na(modDat$TotalTreeCover), "TotalTreeCover"] <- 0.999
# modDat[modDat$CAMCover ==1& !is.na(modDat$CAMCover), "CAMCover"] <- 0.999
# modDat[modDat$TotalHerbaceousCover ==1 & !is.na(modDat$TotalHerbaceousCover), "TotalHerbaceousCover"] <- 0.999
# modDat[modDat$BareGroundCover ==1& !is.na(modDat$BareGroundCover), "BareGroundCover"] <- 0.999
# modDat[modDat$ShrubCover ==1& !is.na(modDat$ShrubCover), "ShrubCover"] <- 0.999
# modDat[modDat$BroadleavedTreeCover_prop ==1& !is.na(modDat$BroadleavedTreeCover_prop), "BroadleavedTreeCover_prop"] <- 0.999
# modDat[modDat$NeedleLeavedTreeCover_prop ==1& !is.na(modDat$NeedleLeavedTreeCover_prop), "NeedleLeavedTreeCover_prop"] <- 0.999
# modDat[modDat$C4Cover_prop ==1& !is.na(modDat$C4Cover_prop), "C4Cover_prop"] <- 0.999
# modDat[modDat$C3Cover_prop ==1& !is.na(modDat$C3Cover_prop), "C3Cover_prop"] <- 0.999
# modDat[modDat$ForbCover_prop ==1& !is.na(modDat$ForbCover_prop), "ForbCover_prop"] <- 0.999

Prep data

set.seed(1234)
modDat_1 <- modDat %>% 
  select(-c(prcp_annTotal:annVPD_min)) %>% 
  # mutate(Lon = st_coordinates(.)[,1], 
  #        Lat = st_coordinates(.)[,2])  %>% 
  # st_drop_geometry() %>% 
  # filter(!is.na(newRegion))
  rename("tmin" = tmin_meanAnnAvg_CLIM, 
     "tmax" = tmax_meanAnnAvg_CLIM, #1
     "tmean" = tmean_meanAnnAvg_CLIM, 
     "prcp" = prcp_meanAnnTotal_CLIM, 
     "t_warm" = T_warmestMonth_meanAnnAvg_CLIM,
     "t_cold" = T_coldestMonth_meanAnnAvg_CLIM, 
     "prcp_wet" = precip_wettestMonth_meanAnnAvg_CLIM,
     "prcp_dry" = precip_driestMonth_meanAnnAvg_CLIM, 
     "prcp_seasonality" = precip_Seasonality_meanAnnAvg_CLIM, #2
     "prcpTempCorr" = PrecipTempCorr_meanAnnAvg_CLIM,  #3
     "abvFreezingMonth" = aboveFreezing_month_meanAnnAvg_CLIM, 
     "isothermality" = isothermality_meanAnnAvg_CLIM, #4
     "annWatDef" = annWaterDeficit_meanAnnAvg_CLIM, 
     "annWetDegDays" = annWetDegDays_meanAnnAvg_CLIM,
     "VPD_mean" = annVPD_mean_meanAnnAvg_CLIM, 
     "VPD_max" = annVPD_max_meanAnnAvg_CLIM, #5
     "VPD_min" = annVPD_min_meanAnnAvg_CLIM, #6
     "VPD_max_95" = annVPD_max_95percentile_CLIM, 
     "annWatDef_95" = annWaterDeficit_95percentile_CLIM, 
     "annWetDegDays_5" = annWetDegDays_5percentile_CLIM, 
     "frostFreeDays_5" = durationFrostFreeDays_5percentile_CLIM, 
     "frostFreeDays" = durationFrostFreeDays_meanAnnAvg_CLIM, 
     "soilDepth" = soilDepth, #7
     "clay" = surfaceClay_perc, 
     "sand" = avgSandPerc_acrossDepth, #8
     "coarse" = avgCoarsePerc_acrossDepth, #9
     "carbon" = avgOrganicCarbonPerc_0_3cm, #10
     "AWHC" = totalAvailableWaterHoldingCapacity,
     ## anomaly variables
     tmean_anom = tmean_meanAnnAvg_3yrAnom, #15
     tmin_anom = tmin_meanAnnAvg_3yrAnom, #16
     tmax_anom = tmax_meanAnnAvg_3yrAnom, #17
    prcp_anom = prcp_meanAnnTotal_3yrAnom, #18
      t_warm_anom = T_warmestMonth_meanAnnAvg_3yrAnom,  #19
     t_cold_anom = T_coldestMonth_meanAnnAvg_3yrAnom, #20
      prcp_wet_anom = precip_wettestMonth_meanAnnAvg_3yrAnom, #21
      precp_dry_anom = precip_driestMonth_meanAnnAvg_3yrAnom,  #22
    prcp_seasonality_anom = precip_Seasonality_meanAnnAvg_3yrAnom, #23 
     prcpTempCorr_anom = PrecipTempCorr_meanAnnAvg_3yrAnom, #24
      aboveFreezingMonth_anom = aboveFreezing_month_meanAnnAvg_3yrAnom, #25  
    isothermality_anom = isothermality_meanAnnAvg_3yrAnom, #26
       annWatDef_anom = annWaterDeficit_meanAnnAvg_3yrAnom, #27
     annWetDegDays_anom = annWetDegDays_meanAnnAvg_3yrAnom,  #28
      VPD_mean_anom = annVPD_mean_meanAnnAvg_3yrAnom, #29
      VPD_min_anom = annVPD_min_meanAnnAvg_3yrAnom,  #30
      VPD_max_anom = annVPD_max_meanAnnAvg_3yrAnom,  #31
     VPD_max_95_anom = annVPD_max_95percentile_3yrAnom, #32
      annWatDef_95_anom = annWaterDeficit_95percentile_3yrAnom, #33 
      annWetDegDays_5_anom = annWetDegDays_5percentile_3yrAnom ,  #34
    frostFreeDays_5_anom = durationFrostFreeDays_5percentile_3yrAnom, #35 
      frostFreeDays_anom = durationFrostFreeDays_meanAnnAvg_3yrAnom #36
  )

# small dataset for if testing the data
if(test_run) {
  modDat_1 <- slice_sample(modDat_1, n = 1e5)
}

Add a constant to the response variable (+1) so that models run…

modDat_1[,response] <- modDat_1[,response]+2

Identify the ecoregion and response variable type to use in this model run

ecoregion <- params$ecoregion
response <- params$response
print(paste0("In this model run, the ecoregion is ", ecoregion," and the response variable is ",response))
## [1] "In this model run, the ecoregion is shrubGrass and the response variable is TotalHerbaceousCover"

Subset the data to only include data for the ecoregion of interest

if (ecoregion == "shrubGrass") {
  # select data for the ecoregion of interest
  modDat_1 <- modDat_1 %>%
    filter(newRegion == "dryShrubGrass")
} else if (ecoregion == "forest") {
  # select data for the ecoregion of interest
  modDat_1 <- modDat_1 %>% 
    filter(newRegion %in% c("eastForest", "westForest"))
}

# remove the rows that have no observations for the response variable of interest
modDat_1 <- modDat_1[!is.na(modDat_1[,response]),]

Currently, subsampling data from the “Texas Coastal Plain”, since it’s quite different from other regions and is really messing with model fit

modDat_1_noLA <- modDat_1 %>% 
  filter(NA_L2NAME != "TEXAS-LOUISIANA COASTAL PLAIN")
modDat_1_LA <- modDat_1 %>% 
  filter(NA_L2NAME == "TEXAS-LOUISIANA COASTAL PLAIN")
# sample points 
modDat_1 <- modDat_1_LA %>% 
  slice_sample(n = round(nrow(modDat_1_LA)*.3)) %>% 
  rbind(modDat_1_noLA) 

Visualize the response variable

hist(modDat_1[,response], main = paste0("Histogram of ",response),
     xlab = paste0(response))

Visualize the predictor variables

The following are the candidate predictor variables for this ecoregion:

if (ecoregion == "shrubGrass") {
  # select potential predictor variables for the ecoregion of interest
        prednames <-
          c(
"tmean"             , "prcp"                    ,"prcp_seasonality"        ,"prcpTempCorr"          , 
"isothermality"     , "annWatDef"               ,"sand"                    ,"coarse"                , 
"carbon"            , "AWHC"                    ,"tmin_anom"               ,"tmax_anom"             , 
"t_warm_anom"       , "prcp_wet_anom"           ,"precp_dry_anom"          ,"prcp_seasonality_anom" , 
"prcpTempCorr_anom" , "aboveFreezingMonth_anom" ,"isothermality_anom"      ,"annWatDef_anom"        , 
"annWetDegDays_anom", "VPD_mean_anom"           ,"VPD_min_anom"            ,"frostFreeDays_5_anom"   )
  
} else if (ecoregion == "forest") {
  # select potential predictor variables for the ecoregion of interest
  prednames <- 
    c(
"tmean"                 ,"prcp"               , "prcp_dry"                , "prcpTempCorr"      ,     
"isothermality"         ,"annWatDef"          , "clay"                    , "sand"              ,     
"coarse"                ,"carbon"             , "AWHC"                    , "tmin_anom"         ,     
"tmax_anom"             ,"prcp_anom"          , "prcp_wet_anom"           , "precp_dry_anom"    ,     
"prcp_seasonality_anom" ,"prcpTempCorr_anom"  , "aboveFreezingMonth_anom" , "isothermality_anom",     
"annWatDef_anom"        ,"annWetDegDays_anom" , "VPD_mean_anom"           , "VPD_max_95_anom"   ,     
"frostFreeDays_5_anom"   )
}

# subset the data to only include these predictors, and remove any remaining NAs 
modDat_1 <- modDat_1 %>% 
  select(prednames, response, newRegion, Year, Long, Lat, NA_L1NAME, NA_L2NAME) %>% 
  drop_na()

names(prednames) <- prednames
df_pred <- modDat_1[, prednames]
# 
# # print the list of predictor variables
# knitr::kable(format = "html", data.frame("Possible_Predictors" = prednames)
# ) %>%
#   kable_styling(bootstrap_options = c("striped", "hover", "condensed")) 
create_summary <- function(df) {
  df %>% 
    pivot_longer(cols = everything(),
                 names_to = 'variable') %>% 
    group_by(variable) %>% 
    summarise(across(value, .fns = list(mean = ~mean(.x, na.rm = TRUE), min = ~min(.x, na.rm = TRUE), 
                                        median = ~median(.x, na.rm = TRUE), max = ~max(.x, na.rm = TRUE)))) %>% 
    mutate(across(where(is.numeric), round, 4))
}

modDat_1[prednames] %>% 
  create_summary() %>% 
  knitr::kable(caption = 'summaries of possible predictor variables') %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed")) 
summaries of possible predictor variables
variable value_mean value_min value_median value_max
AWHC 14.4020 0.6319 13.8921 35.1058
VPD_mean_anom -0.0200 -0.3600 -0.0190 0.2199
VPD_min_anom -0.0125 -0.1764 -0.0158 0.1390
aboveFreezingMonth_anom 0.0556 -1.3667 0.0370 2.2143
annWatDef 99.3328 1.5871 88.0773 451.1875
annWatDef_anom -0.0636 -3.9347 -0.0617 1.0000
annWetDegDays_anom 0.0070 -1.6606 0.0074 0.9840
carbon 1.0617 0.0001 0.8324 34.1410
coarse 8.2852 0.0000 5.3244 77.8683
frostFreeDays_5_anom -16.7387 -244.0000 -6.0000 55.8000
isothermality 37.4480 21.6016 37.2058 58.2327
isothermality_anom 0.4811 -6.0230 0.4170 8.9119
prcp 420.9280 90.8317 372.6317 1690.7422
prcpTempCorr 0.0556 -0.7786 0.1111 0.6890
prcpTempCorr_anom 0.0168 -0.6495 0.0195 0.5744
prcp_seasonality 0.9356 0.5255 0.9205 1.7468
prcp_seasonality_anom -0.0162 -0.8313 -0.0082 0.4743
prcp_wet_anom -0.0112 -1.9576 0.0087 0.6958
precp_dry_anom -0.0409 -9.0000 0.2508 1.0000
sand 46.0031 0.5324 45.6216 98.9991
t_warm_anom -0.4626 -6.3974 -0.4711 3.3598
tmax_anom -0.2707 -3.4286 -0.3059 2.3198
tmean 10.5471 -2.2935 9.3595 23.9289
tmin_anom -0.4508 -4.0550 -0.4337 2.4361
# response_summary <- modDat_1 %>% 
#     dplyr::select(#where(is.numeric), -all_of(pred_vars),
#       matches(response)) %>% 
#     create_summary()
# 
# 
# kable(response_summary, 
#       caption = 'summaries of response variables, calculated using paint') %>%
# kable_styling(bootstrap_options = c("striped", "hover", "condensed")) 

Plot predictor vars against each other

set.seed(12011993)
# function for colors
my_fn <- function(data, mapping, method="p", use="pairwise", ...){
  
  # grab data
  x <- eval_data_col(data, mapping$x)
  y <- eval_data_col(data, mapping$y)
  
  # calculate correlation
  corr <- cor(x, y, method=method, use=use)
  
  # calculate colour based on correlation value
  # Here I have set a correlation of minus one to blue, 
  # zero to white, and one to red 
  # Change this to suit: possibly extend to add as an argument of `my_fn`
  colFn <- colorRampPalette(c("red", "white", "blue"), interpolate ='spline')
  fill <- colFn(100)[findInterval(corr, seq(-1, 1, length=100))]
  
  ggally_cor(data = data, mapping = mapping, size = 2.5, stars = FALSE, 
             digits = 2, colour = I("black"),...) + 
    theme_void() +
    theme(panel.background = element_rect(fill=fill))
  
}

if (run == TRUE) {
(corrPlot <- modDat_1 %>% 
  select(prednames) %>% 
  slice_sample(n = 5e4) %>% 
  #select(-matches("_")) %>% 
ggpairs( upper = list(continuous = my_fn, size = .1), lower = list(continuous = GGally::wrap("points", alpha = 0.1, size=0.1)), progress = FALSE))
    base::saveRDS(corrPlot, paste0("../ModelFitting/models/", response, "_",ecoregion, "_corrPlot.rds"))
  
  } else {
    # corrPlot <- readRDS(paste0("../ModelFitting/models/", response, "_",ecoregion, "_corrPlot.rds"))
    # (corrPlot)
    print(c("See previous correlation figures"))
  }

Predictor variables compared to binned response variables

set.seed(12011993)
# vector of name of response variables
vars_response <- response

# longformat dataframes for making boxplots
df_sample_plots <-  modDat_1  %>% 
  slice_sample(n = 5e4) %>% 
   rename(response = all_of(response)) %>% 
  mutate(response = case_when(
    response <= .25 ~ ".25", 
    response > .25 & response <=.5 ~ ".5", 
    response > .5 & response <=.75 ~ ".75", 
    response >= .75  ~ "1", 
  )) %>% 
  select(c(response, prednames)) %>% 
  tidyr::pivot_longer(cols = unname(prednames), 
               names_to = "predictor", 
               values_to = "value"
               )  
 

  ggplot(df_sample_plots, aes_string(x= "response", y = 'value')) +
  geom_boxplot() +
  facet_wrap(~predictor , scales = 'free_y') + 
  ylab("Predictor Variable Values") + 
    xlab(response)

Standardize the predictor variables for the model-fitting process

modDat_1_s <- modDat_1 %>% 
  mutate(across(all_of(prednames), base::scale, .names = "{.col}_s")) 
names(modDat_1_s) <- c(names(modDat_1),
                       paste0(prednames, "_s")
                       )
  
scaleFigDat_1 <- modDat_1_s %>% 
  select(c(Long, Lat, Year, prednames)) %>% 
  pivot_longer(cols = all_of(names(prednames)), 
               names_to = "predNames", 
               values_to = "predValues_unScaled")
scaleFigDat_2 <- modDat_1_s %>% 
  select(c(Long, Lat, Year,paste0(prednames, "_s"))) %>% 
  pivot_longer(cols = all_of(paste0(prednames, "_s")), 
               names_to = "predNames", 
               values_to = "predValues_scaled", 
               names_sep = ) %>% 
  mutate(predNames = str_replace(predNames, pattern = "_s$", replacement = ""))

scaleFigDat_3 <- scaleFigDat_1 %>% 
  left_join(scaleFigDat_2)

ggplot(scaleFigDat_3) + 
  facet_wrap(~predNames, scales = "free") +
  geom_histogram(aes(predValues_unScaled), fill = "lightgrey", col = "darkgrey") + 
  geom_histogram(aes(predValues_scaled), fill = "lightblue", col = "blue") +
  xlab ("predictor variable values") + 
  ggtitle("Comparing the distribution of unscaled (grey) to scaled (blue) predictor variables")

Model Fitting

Visualize the level 2 ecoregions and how they differ across environmental space

## visualize the variation between groups across environmental space

## do a pca of climate across level 2 ecoregions
pca <- prcomp(modDat_1_s[,paste0(prednames, "_s")])
library(factoextra)
(fviz_pca_ind(pca, habillage = modDat_1_s$NA_L2NAME, label = "none", addEllipses = TRUE, ellipse.level = .95, ggtheme = theme_minimal(), alpha.ind = .1))

if (ecoregion == "shrubGrass") {
  print("We'll combine the 'Mediterranean California' and 'Western Sierra Madre Piedmont' ecoregions (into 'Mediterranean Piedmont'). We'll also combine `Tamaulipas-Texas semiarid plain,' 'Texas-Lousiana Coastal plain,' and 'South Central semiarid prairies' ecoregions (into (`Semiarid plain and prairies`)." )
  
  modDat_1_s[modDat_1_s$NA_L2NAME %in% c("MEDITERRANEAN CALIFORNIA", "WESTERN SIERRA MADRE PIEDMONT"), "NA_L2NAME"] <- "MEDITERRANEAN PIEDMONT"
  modDat_1[modDat_1$NA_L2NAME %in% c("MEDITERRANEAN CALIFORNIA", "WESTERN SIERRA MADRE PIEDMONT"), "NA_L2NAME"] <- "MEDITERRANEAN PIEDMONT"
  
  modDat_1_s[modDat_1_s$NA_L2NAME %in% c("TAMAULIPAS-TEXAS SEMIARID PLAIN", "TEXAS-LOUISIANA COASTAL PLAIN", "SOUTH CENTRAL SEMIARID PRAIRIES"), "NA_L2NAME"] <- "SEMIARID PLAIN AND PRAIRIES"
  modDat_1[modDat_1$NA_L2NAME %in% c("TAMAULIPAS-TEXAS SEMIARID PLAIN", "TEXAS-LOUISIANA COASTAL PLAIN", "SOUTH CENTRAL SEMIARID PRAIRIES"), "NA_L2NAME"] <- "SEMIARID PLAIN AND PRAIRIES"
}
## [1] "We'll combine the 'Mediterranean California' and 'Western Sierra Madre Piedmont' ecoregions (into 'Mediterranean Piedmont'). We'll also combine `Tamaulipas-Texas semiarid plain,' 'Texas-Lousiana Coastal plain,' and 'South Central semiarid prairies' ecoregions (into (`Semiarid plain and prairies`)."
# make a table of n for each region
modDat_1 %>% 
  group_by(NA_L2NAME) %>% 
  dplyr::summarize("Number_Of_Observations" = length(NA_L2NAME)) %>% 
  rename("Level_2_Ecoregion" = NA_L2NAME)%>% 
  kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed")) 
Level_2_Ecoregion Number_Of_Observations
COLD DESERTS 64212
MEDITERRANEAN PIEDMONT 2897
SEMIARID PLAIN AND PRAIRIES 25614
TEMPERATE PRAIRIES 3314
WARM DESERTS 10186
WEST-CENTRAL SEMIARID PRAIRIES 19694

Then, look at the spatial distribution and environmental characteristics of the grouped ecoregions

## make data into spatial format
modDat_1_sf <- modDat_1 %>% 
  st_as_sf(coords = c("Long", "Lat"), crs = st_crs("PROJCRS[\"unnamed\",\n    BASEGEOGCRS[\"unknown\",\n        DATUM[\"unknown\",\n            ELLIPSOID[\"Spheroid\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1,\n                    ID[\"EPSG\",9001]]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433,\n                ID[\"EPSG\",9122]]]],\n    CONVERSION[\"Lambert Conic Conformal (2SP)\",\n        METHOD[\"Lambert Conic Conformal (2SP)\",\n            ID[\"EPSG\",9802]],\n        PARAMETER[\"Latitude of false origin\",42.5,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8821]],\n        PARAMETER[\"Longitude of false origin\",-100,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8822]],\n        PARAMETER[\"Latitude of 1st standard parallel\",25,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8823]],\n        PARAMETER[\"Latitude of 2nd standard parallel\",60,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8824]],\n        PARAMETER[\"Easting at false origin\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8826]],\n        PARAMETER[\"Northing at false origin\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8827]]],\n    CS[Cartesian,2],\n        AXIS[\"easting\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]],\n        AXIS[\"northing\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]]]"))


# download map info for visualization
data(state_boundaries_wgs84) 

cropped_states <- suppressMessages(state_boundaries_wgs84 %>%
  dplyr::filter(NAME!="Hawaii") %>%
  dplyr::filter(NAME!="Alaska") %>%
  dplyr::filter(NAME!="Puerto Rico") %>%
  dplyr::filter(NAME!="American Samoa") %>%
  dplyr::filter(NAME!="Guam") %>%
  dplyr::filter(NAME!="Commonwealth of the Northern Mariana Islands") %>%
  dplyr::filter(NAME!="United States Virgin Islands") %>%

  sf::st_sf() %>%
  sf::st_transform(sf::st_crs(modDat_1_sf))) #%>%
  #sf::st_crop(sf::st_bbox(modDat_1_sf)+c(-1,-1,1,1))

map1 <- ggplot() +
  geom_sf(data=cropped_states,fill='white') +
  geom_sf(data=modDat_1_sf,aes(fill=as.factor(NA_L2NAME)),linewidth=0.5,alpha=0.5) +
  geom_point(data=modDat_1,alpha=0.5, 
             aes(x = Long, y = Lat, color=as.factor(NA_L2NAME)), alpha = .1) +
  #scale_fill_okabeito() +
  #scale_color_okabeito() +
 # theme_default() +
  theme(legend.position = 'none') +
  labs(title = "Level 2 Ecoregions as spatial blocks")

hull <- modDat_1_sf %>%
  ungroup() %>%
  group_by(NA_L2NAME) %>%
  slice(chull(tmean, prcp))

plot1<-ggplot(data=modDat_1_sf,aes(x=tmean,y=prcp)) +
  geom_polygon(data = hull, alpha = 0.25,aes(fill=NA_L2NAME) )+
  geom_point(aes(group=NA_L2NAME,color=NA_L2NAME),alpha=0.25) +
  theme_minimal() + xlab("Annual Average T_mean - long-term average") +
  ylab("Annual Average Precip - long-term average") #+
  #scale_color_okabeito() +
  #scale_fill_okabeito()

plot2<-ggplot(data=modDat_1_sf %>%
                pivot_longer(cols=tmean:prcp),
              aes(x=value,group=name)) +
  # geom_polygon(data = hull, alpha = 0.25,aes(fill=fold) )+
  geom_density(aes(group=NA_L2NAME,fill=NA_L2NAME),alpha=0.25) +
  theme_minimal() +
  facet_wrap(~name,scales='free')# +
  #scale_color_okabeito() +
  #scale_fill_okabeito()
 
library(patchwork)
(combo <- (map1+plot1)/plot2) 

Fit a global model with all of the data

First, fit a LASSO regression model using the glmnet R package

  • This regression is a Gamma glm with a log link
  • Use cross validation across level 2 ecoregions to tune the lambda parameter in the LASSO model
  • this model is fit to using the scaled weather/climate/soils variables
  • this list of possible predictors includes:
    1. main effects
    2. interactions between all soils variables
    3. interactions between climate and weather variables
    4. transformed main effects (squared, log-transformed (add a uniform integer – 20– to all variables prior to log-transformation), square root-transformed (add a uniform integer – 20– to all variables prior to log-transformation))
## 
## Call:  cv.glmnet(x = X[, 2:ncol(X)], y = y, type.measure = "mse", foldid = my_folds,      keep = TRUE, parallel = TRUE, family = stats::Gamma(link = "log"),      alpha = 1, nlambda = 100, standardize = FALSE) 
## 
## Measure: Mean-Squared Error 
## 
##     Lambda Index Measure    SE Nonzero
## min  0.011    39    1076 366.2      64
## 1se  0.379     1    1230 145.5       0

Then, fit regular glm models (Gamma glm with a log link), first using the coefficients from the ‘best’ lambda identified in the LASSO model, as then using the coefficients from the ‘1SE’ lambda identified from the LASSO (this is the value of lambda such that the cross validation error is within 1 standard error of the minimum).

## fit w/ the identified coefficients from the 'best' lambda, but using the glm function
  mat_glmnet_best <- as.matrix(bestLambda_coef)
  mat2_glmnet_best <- mat_glmnet_best[mat_glmnet_best[,1] != 0,]
  names(mat2_glmnet_best) <- rownames(mat_glmnet_best)[mat_glmnet_best[,1] != 0]

  if (length(mat2_glmnet_best) == 1) {
    f_glm_bestLambda <- as.formula(paste0(response, "~ 1"))
  } else {
  f_glm_bestLambda <- as.formula(paste0(response, " ~ ", paste0(names( mat2_glmnet_best)[2:length(names( mat2_glmnet_best))], collapse = " + ")))
  }
  
  fit_glm_bestLambda <- glm(data = modDat_1_s
                              , formula =  f_glm_bestLambda, family =  stats::Gamma(link = "log"))
  
   ## fit w/ the identified coefficients from the '1se' lambda, but using the glm function
  mat_glmnet_1se <- as.matrix(seLambda_coef)
  mat2_glmnet_1se <- mat_glmnet_1se[mat_glmnet_1se[,1] != 0,]
  names(mat2_glmnet_1se) <- rownames(mat_glmnet_1se)[mat_glmnet_1se[,1] != 0]
  if(length(mat2_glmnet_1se) == 1) {
    f_glm_1se <- as.formula(paste0(response, "~ 1"))
  } else {
  f_glm_1se <- as.formula(paste0(response, " ~ ", paste0(names( mat2_glmnet_1se)[2:length(names( mat2_glmnet_1se))], collapse = " + ")))
  }


  fit_glm_se <- glm(data = modDat_1_s, formula = f_glm_1se,
                    family =  stats::Gamma(link = "log"))

Then, we predict (on the training set) using both of these models (best lambda and 1se lambda)

  ## predict on the test data
  # lasso model predictions with the optimal lambda
  optimal_pred <- predict(fit_glm_bestLambda, newx=X[,2:ncol(X)], type = "response") - 2
  optimal_pred_1se <-  predict(fit_glm_se, newx=X[,2:ncol(X)], type = "response") - 2
    null_fit <- glm(#data = data.frame("y" = y, X[,paste0(prednames, "_s")]), 
      formula = y ~ 1, family = stats::Gamma(link = "log"))
  null_pred <- predict(null_fit, newdata = as.data.frame(X), type = "response"
                       ) - 2

  # save data
  fullModOut <- list(
    "modelObject" = fit,
    "nullModelObject" = null_fit,
    "modelPredictions" = data.frame(#ecoRegion_holdout = rep(test_eco,length(y)),
      obs=y,
                    pred_opt=optimal_pred, 
                    pred_opt_se = optimal_pred_1se,
                    pred_null=null_pred#,
                    #pred_nopenalty=nopen_pred
                    ))
  
  
# calculate correlations between null and optimal model 
my_cors <- c(cor(optimal_pred, c(y-2)),
             cor(optimal_pred_1se, c(y-2)), 
            cor(null_pred, c(y-2))
            )

# calculate mse between null and optimal model 
my_mse <- c(mean((fullModOut$modelPredictions$pred_opt -  c(y-2))^2) ,
            mean((fullModOut$modelPredictions$pred_opt_se -  c(y-2))^2) ,
            mean((fullModOut$modelPredictions$pred_null - c(y-2))^2)#,
            #mean((obs_pred$pred_nopenalty - obs_pred$obs)^2)
            )

ggplot() + 
  geom_point(aes(X[,2], fullModOut$modelPredictions$obs-2), col = "black", alpha = .1) + 
  geom_point(aes(X[,2], fullModOut$modelPredictions$pred_opt), col = "red", alpha = .1) + ## predictions w/ the CV model
  geom_point(aes(X[,2], fullModOut$modelPredictions$pred_opt_se), col = "green", alpha = .1) + ## predictions w/ the CV model (1se lambda)
  geom_point(aes(X[,2], fullModOut$modelPredictions$pred_null), col = "blue", alpha = .1) + 
  labs(title = "A rough comparison of observed and model-predicted values", 
       subtitle = "black = observed values \n red = predictions from 'best lambda' model \n green = predictions from '1se' lambda model \n blue = predictions from null model") +
  xlab(colnames(X)[2])

  #ylim(c(0,200))

The internal cross-validation process to fit the global LASSO model identified an optimal lambda value (regularization parameter) of r{print(best_lambda)}. The lambda value such that the cross validation error is within 1 standard error of the minimum (“1se lambda”) was `r{print(fit$lambda.1se)}`` . The following coefficients were kept in each model:

# the coefficient matrix from the 'best model' -- find and print those coefficients that aren't 0 in a table
coef_glm_bestLambda <- coef(fit_glm_bestLambda) %>% 
  data.frame() 
coef_glm_bestLambda$coefficientName <- rownames(coef_glm_bestLambda)
names(coef_glm_bestLambda)[1] <- "coefficientValue_bestLambda"
# coefficient matrix from the '1se' model 
coef_glm_1se <- coef(fit_glm_se) %>% 
  data.frame() 
coef_glm_1se$coefficientName <- rownames(coef_glm_1se)
names(coef_glm_1se)[1] <- "coefficientValue_1seLambda"
# add together
coefs <- full_join(coef_glm_bestLambda, coef_glm_1se) %>% 
  select(coefficientName, coefficientValue_bestLambda, coefficientValue_1seLambda)

globModTerms <- coefs[!is.na(coefs$coefficientValue_bestLambda), "coefficientName"]

## also, get the number of unique variables in each model 
var_prop_pred <- paste0(response, "_pred")
response_vars <- c(response, var_prop_pred)
# for best lambda model
prednames_fig <- paste(str_split(globModTerms, ":", simplify = TRUE)) 
prednames_fig <- str_replace(prednames_fig, "I\\(", "")
prednames_fig <- str_replace(prednames_fig, "\\^2\\)", "")
prednames_fig <- unique(prednames_fig[prednames_fig>0])
prednames_fig <- prednames_fig
prednames_fig_num <- length(prednames_fig)
# for 1SE lambda model
globModTerms_1se <- coefs[!is.na(coefs$coefficientValue_1seLambda), "coefficientName"]
if (length(globModTerms_1se) == 1) {
prednames_fig_1se <- paste(str_split(globModTerms_1se, ":", simplify = TRUE)) 
prednames_fig_1se <- str_replace(prednames_fig_1se, "I\\(", "")
prednames_fig_1se <- str_replace(prednames_fig_1se, "\\^2\\)", "")
prednames_fig_1se <- unique(prednames_fig_1se[prednames_fig_1se>0])
prednames_fig_1se_num <- c(0)
} else {
prednames_fig_1se <- paste(str_split(globModTerms_1se, ":", simplify = TRUE)) 
prednames_fig_1se <- str_replace(prednames_fig_1se, "I\\(", "")
prednames_fig_1se <- str_replace(prednames_fig_1se, "\\^2\\)", "")
prednames_fig_1se <- unique(prednames_fig_1se[prednames_fig_1se>0])
prednames_fig_1se_num <- length(prednames_fig_1se)
}

# make a table
kable(coefs, col.names = c("Coefficient Name", "Value from best lambda model", "Value from 1se lambda model")
      ) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed")) 
Coefficient Name Value from best lambda model Value from 1se lambda model
(Intercept) 3.4100497 3.660819
prcp_s 0.2106097 NA
prcpTempCorr_s 0.2592503 NA
isothermality_s -0.1176705 NA
annWatDef_s -0.0088556 NA
sand_s -0.0913288 NA
coarse_s -0.0133341 NA
AWHC_s 0.1075600 NA
t_warm_anom_s 0.0177688 NA
isothermality_anom_s -0.0114444 NA
annWetDegDays_anom_s -0.0137544 NA
I(prcp_s^2) -0.0780247 NA
I(prcpTempCorr_s^2) 0.1862365 NA
I(isothermality_s^2) -0.0282653 NA
I(tmin_anom_s^2) -0.0212262 NA
I(prcp_seasonality_anom_s^2) -0.0038217 NA
I(prcpTempCorr_anom_s^2) -0.0125591 NA
I(aboveFreezingMonth_anom_s^2) 0.0103105 NA
I(annWetDegDays_anom_s^2) -0.0220411 NA
I(VPD_mean_anom_s^2) -0.0120781 NA
I(VPD_min_anom_s^2) 0.0211411 NA
I(sand_s^2) 0.0293809 NA
I(carbon_s^2) 0.0018767 NA
I(AWHC_s^2) -0.0422001 NA
annWetDegDays_anom_s:aboveFreezingMonth_anom_s -0.0175734 NA
aboveFreezingMonth_anom_s:prcp_wet_anom_s -0.0059107 NA
aboveFreezingMonth_anom_s:prcpTempCorr_anom_s 0.0242005 NA
aboveFreezingMonth_anom_s:VPD_min_anom_s 0.0162988 NA
prcp_s:annWatDef_anom_s -0.0085184 NA
annWatDef_s:frostFreeDays_5_anom_s -0.0111776 NA
prcpTempCorr_s:annWatDef_s -0.1719310 NA
annWatDef_s:tmin_anom_s -0.0298934 NA
annWetDegDays_anom_s:frostFreeDays_5_anom_s 0.0080789 NA
annWetDegDays_anom_s:prcpTempCorr_anom_s -0.0059790 NA
prcpTempCorr_s:annWetDegDays_anom_s 0.0147086 NA
annWetDegDays_anom_s:tmax_anom_s 0.0026184 NA
isothermality_anom_s:frostFreeDays_5_anom_s 0.0125496 NA
prcp_wet_anom_s:frostFreeDays_5_anom_s 0.0052059 NA
frostFreeDays_5_anom_s:tmean_s -0.0152895 NA
isothermality_anom_s:tmin_anom_s -0.0156470 NA
isothermality_anom_s:VPD_min_anom_s 0.0043528 NA
prcp_s:isothermality_s 0.0696158 NA
isothermality_s:prcp_seasonality_s 0.0704052 NA
prcpTempCorr_s:isothermality_s 0.0688654 NA
isothermality_s:t_warm_anom_s 0.0041256 NA
isothermality_s:tmax_anom_s 0.0225992 NA
isothermality_s:tmin_anom_s 0.0378574 NA
prcp_s:tmean_s 0.0876542 NA
prcp_seasonality_s:prcp_seasonality_anom_s 0.0178788 NA
prcpTempCorr_anom_s:prcp_seasonality_s -0.0120024 NA
prcpTempCorr_s:prcp_seasonality_s 0.0842187 NA
t_warm_anom_s:prcp_wet_anom_s -0.0309288 NA
prcp_wet_anom_s:tmax_anom_s 0.0327531 NA
prcp_wet_anom_s:tmean_s 0.0094295 NA
prcp_wet_anom_s:VPD_min_anom_s -0.0270501 NA
t_warm_anom_s:prcpTempCorr_anom_s 0.0177379 NA
prcpTempCorr_s:precp_dry_anom_s 0.0071519 NA
prcpTempCorr_s:tmax_anom_s 0.0248815 NA
prcpTempCorr_s:tmin_anom_s 0.0108157 NA
tmean_s:precp_dry_anom_s 0.0120453 NA
tmin_anom_s:tmax_anom_s 0.0297353 NA
VPD_min_anom_s:tmean_s -0.0083570 NA
coarse_s:AWHC_s 0.0195495 NA
sand_s:AWHC_s 0.0561811 NA
sand_s:coarse_s -0.0157758 NA
# calculate RMSE of both models 
RMSE_best <- yardstick::rmse(fullModOut$modelPredictions[,c("obs", "pred_opt")]-2, truth = "obs", estimate = "pred_opt")$.estimate
RMSE_1se <- yardstick::rmse(fullModOut$modelPredictions[,c("obs", "pred_opt_se")]-2, truth = "obs", estimate = "pred_opt_se")$.estimate
bias_best <-  mean(fullModOut$modelPredictions$obs - fullModOut$modelPredictions$pred_opt)
bias_1se <- mean(fullModOut$modelPredictions$obs - fullModOut$modelPredictions$pred_opt_se)

uniqueCoeffs <- data.frame("Best lambda model" = c(RMSE_best, bias_best,
  as.integer(length(globModTerms)-1), as.integer(prednames_fig_num), 
                                                   as.integer(sum(prednames_fig %in% c(prednames_clim))),
                                                   as.integer(sum(prednames_fig %in% c(prednames_weath))),
                                                   as.integer(sum(prednames_fig %in% c(prednames_soils)))
                                                   ), 
                           "1se lambda model" = c(RMSE_1se, bias_1se,
                             length(globModTerms_1se)-1, prednames_fig_1se_num,
                                                   sum(prednames_fig_1se %in% c(prednames_clim)),
                                                   sum(prednames_fig_1se %in% c(prednames_weath)),
                                                   sum(prednames_fig_1se %in% c(prednames_soils))))
row.names(uniqueCoeffs) <- c("RMSE", "bias - mean(obs-pred.)", "Total number of coefficients", "Number of unique coefficients",
                             "Number of unique climate coefficients", 
                             "Number of unique weather coefficients",  
                             "Number of unique soils coefficients"
                             )

kable(uniqueCoeffs, 
      col.names = c("Best lambda model", "1se lambda model"), row.names = TRUE) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed")) 
Best lambda model 1se lambda model
RMSE 20.880466 28.8795
bias - mean(obs-pred.) 1.846602 2.0000
Total number of coefficients 64.000000 0.0000
Number of unique coefficients 24.000000 0.0000
Number of unique climate coefficients 6.000000 0.0000
Number of unique weather coefficients 14.000000 0.0000
Number of unique soils coefficients 4.000000 0.0000

Visualizations of Model Predictions and Residuals – using best lambda model

observed vs. predicted values

Predicting on the data

  # create prediction for each each model
# (i.e. for each fire proporation variable)
predict_by_response <- function(mod, df) {
  df_out <- df
  response_name <- paste0(response, "_pred")
  df_out <- df_out %>% cbind(predict(mod, newx= df_out, #s="lambda.min", 
                                     type = "response"))
   colnames(df_out)[ncol(df_out)] <- response_name
  return(df_out)
}

pred_glm1 <- predict_by_response(fit_glm_bestLambda, X[,2:ncol(X)])

# add back in true y values
pred_glm1 <- pred_glm1 %>% 
  cbind( data.frame("y" = y))
# rename the true response column to not be 'y_test' 
colnames(pred_glm1)[which(colnames(pred_glm1) == "y")] <- paste(response)

# add back in lat/long data 
pred_glm1 <- pred_glm1 %>% 
  cbind(modDat_1_s[,c("Long", "Lat", "Year")])

pred_glm1$resid <- pred_glm1[,response] - pred_glm1[,paste0(response, "_pred")]
pred_glm1$extremeResid <- NA
pred_glm1[pred_glm1$resid > 70 | pred_glm1$resid < -70,"extremeResid"] <- 1

# plot(x = pred_glm1[,response],
#      y = pred_glm1[,paste0(response, "_pred")],
#      xlab = "observed values", ylab = "predicted values")
# points(x = pred_glm1[!is.na(pred_glm1$extremeResid), response],
#        y = pred_glm1[!is.na(pred_glm1$extremeResid), paste0(response, "_pred")],
#        col = "red")
pred_glm1_1se <- predict_by_response(fit_glm_se, X[,2:ncol(X)])

# add back in true y values
pred_glm1_1se <- pred_glm1_1se %>% 
  cbind( data.frame("y" = y))
# rename the true response column to not be 'y_test' 
colnames(pred_glm1_1se)[which(colnames(pred_glm1_1se) == "y")] <- paste(response)

# add back in lat/long data 
pred_glm1_1se <- pred_glm1_1se %>% 
  cbind(modDat_1_s[,c("Long", "Lat", "Year")])

pred_glm1_1se$resid <- pred_glm1_1se[,response] - pred_glm1_1se[,paste0(response, "_pred")]
pred_glm1_1se$extremeResid <- NA
pred_glm1_1se[pred_glm1_1se$resid > 70 | pred_glm1_1se$resid < -70,"extremeResid"] <- 1

Maps of Observations, Predictions, and Residuals=

Observations across the temporal range of the dataset

pred_glm1 <- pred_glm1 %>% 
  mutate(resid = .[[response]] - .[[paste0(response,"_pred")]]) 

# rasterize
# get reference raster
test_rast <-  rast("../../../Data_raw/dayMet/rawMonthlyData/orders/70e0da02b9d2d6e8faa8c97d211f3546/Daymet_Monthly_V4R1/data/daymet_v4_prcp_monttl_na_1980.tif") %>% 
  terra::aggregate(fact = 8, fun = "mean")
## |---------|---------|---------|---------|=========================================                                          
## add ecoregion boundaries (for our ecoregion level model)
regions <- sf::st_read(dsn = "../../../Data_raw/Level2Ecoregions/", layer = "NA_CEC_Eco_Level2") 
## Reading layer `NA_CEC_Eco_Level2' from data source 
##   `/Users/astears/Documents/Dropbox_static/Work/NAU_USGS_postdoc/PED_vegClimModels/Data_raw/Level2Ecoregions' using driver `ESRI Shapefile'
## Simple feature collection with 2261 features and 8 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -4334052 ymin: -3313739 xmax: 3324076 ymax: 4267265
## Projected CRS: Sphere_ARC_INFO_Lambert_Azimuthal_Equal_Area
regions <- regions %>% 
  st_transform(crs = st_crs(test_rast)) %>% 
  st_make_valid() #%>% 
  #st_crop(st_bbox(test_rast))
# 
# goodRegions_temp <- st_overlaps(y = cropped_states, x = regions, sparse = FALSE) %>% 
#   rowSums() 
# goodRegions <- regions[goodRegions_temp != 0,]

ecoregionLU <- data.frame("NA_L1NAME" = sort(unique(regions$NA_L1NAME)), 
                        "newRegion" = c(NA, "Forest", "dryShrubGrass", 
                                        "dryShrubGrass", "Forest", "dryShrubGrass",
                                       "dryShrubGrass", "Forest", "Forest", 
                                       "dryShrubGrass", "Forest", "Forest", 
                                       "Forest", "Forest", "dryShrubGrass", 
                                       NA
                                        ))
goodRegions <- regions %>% 
  left_join(ecoregionLU)
mapRegions <- goodRegions %>% 
  filter(!is.na(newRegion)) %>% 
  group_by(newRegion) %>% 
  summarise(geometry = sf::st_union(geometry)) %>% 
  ungroup() %>% 
  st_simplify(dTolerance = 1000)
#mapview(mapRegions)
# rasterize data
plotObs <- pred_glm1 %>% 
         drop_na(paste(response)) %>% 
  #slice_sample(n = 5e4) %>%
  terra::vect(geom = c("Long", "Lat")) %>% 
  terra::set.crs(crs(test_rast)) %>% 
  terra::rasterize(y = test_rast, 
                   field = response, 
                   fun = mean) #%>% 
  #terra::aggregate(fact = 2, fun = mean, na.rm = TRUE) %>% 
  #terra::crop(ext(-1950000, 1000000, -1800000, 1000000))

# get the extent of this particular raster, and crop it accordingly
tempExt <- crds(plotObs, na.rm = TRUE)

plotObs_2 <- plotObs %>% 
  crop(ext(min(tempExt[,1]), max(tempExt[,1]),
           min(tempExt[,2]), max(tempExt[,2])) 
       )
# make figures
ggplot() +
geom_spatraster(data = plotObs_2) + 
  geom_sf(data = mapRegions, fill = NA, col = "rosybrown4", lwd = .5) +
  geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = paste0("Observations of ", response, " in the ",ecoregion, " ecoregion")) +
  scale_fill_gradient2(low = "brown",
                       mid = "wheat" ,
                       high = "darkgreen" , 
                       midpoint = 0,   na.value = "lightgrey") + 
  xlim(st_bbox(plotObs_2)[c(1,3)]) + 
  ylim(st_bbox(plotObs_2)[c(2,4)])

Predictions across the temporal range of the dataset

# rasterize data
plotPred <- pred_glm1 %>% 
         drop_na(paste0(response,"_pred")) %>% 
  #slice_sample(n = 5e4) %>%
  terra::vect(geom = c("Long", "Lat")) %>% 
  terra::set.crs(crs(test_rast)) %>% 
  terra::rasterize(y = test_rast, 
                   field = paste0(response,"_pred"), 
                   fun = mean) #%>% 
  #terra::aggregate(fact = 2, fun = mean, na.rm = TRUE) %>% 
  #terra::crop(ext(-1950000, 1000000, -1800000, 1000000))

# get the point location of those predictions that are > 100
highPred_points <- pred_glm1 %>% 
  filter(.[[paste0(response,"_pred")]] > 100 | 
           .[[paste0(response, "_pred")]] < 0) %>% 
  terra::vect(geom = c("Long", "Lat")) %>% 
  terra::set.crs(crs(test_rast)) 

# get the extent of this particular raster, and crop it accordingly
tempExt <- crds(plotPred, na.rm = TRUE)

plotPred_2 <- plotPred %>% 
  crop(ext(min(tempExt[,1]), max(tempExt[,1]),
           min(tempExt[,2]), max(tempExt[,2])) 
       )
# make figures
ggplot() +
geom_spatraster(data = plotPred_2) + 
  geom_sf(data = mapRegions, fill = NA, col = "rosybrown4", lwd = .5) +
  geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_2)),fill=NA )  + 
  geom_sf(data = highPred_points, col = "red") +
labs(title = paste0("Predictions from the 'best lambda' fitted model of ", response, " in the ",ecoregion, " ecoregion"),
     subtitle =  "bestLambda model")  +
  scale_fill_gradient2(low = "wheat",
                       mid = "darkgreen",
                       high = "red" , 
                       midpoint = 100,   na.value = "lightgrey",
                       limits = c(0,100)) + 
  xlim(st_bbox(plotObs_2)[c(1,3)]) + 
  ylim(st_bbox(plotObs_2)[c(2,4)])

# rasterize data
plotPred <- pred_glm1_1se %>% 
         drop_na(paste0(response,"_pred")) %>% 
  #slice_sample(n = 5e4) %>%
  terra::vect(geom = c("Long", "Lat")) %>% 
  terra::set.crs(crs(test_rast)) %>% 
  terra::rasterize(y = test_rast, 
                   field = paste0(response,"_pred"), 
                   fun = mean) #%>% 
  #terra::aggregate(fact = 2, fun = mean, na.rm = TRUE) %>% 
  #terra::crop(ext(-1950000, 1000000, -1800000, 1000000))

# get the point location of those predictions that are > 100
highPred_points <- pred_glm1_1se %>% 
  filter(.[[paste0(response,"_pred")]] > 100 | 
           .[[paste0(response, "_pred")]] < 0) %>% 
  terra::vect(geom = c("Long", "Lat")) %>% 
  terra::set.crs(crs(test_rast)) 

# get the extent of this particular raster, and crop it accordingly
tempExt <- crds(plotPred, na.rm = TRUE)

plotPred_2 <- plotPred %>% 
  crop(ext(min(tempExt[,1]), max(tempExt[,1]),
           min(tempExt[,2]), max(tempExt[,2])) 
       )
# make figures
ggplot() +
geom_spatraster(data = plotPred_2) + 
  geom_sf(data = mapRegions, fill = NA, col = "rosybrown4", lwd = .5) +
  geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_2)),fill=NA )  + geom_sf(data = highPred_points, col = "red") +
labs(title = paste0("Predictions from the '1SE lambda' fitted model of ", response, " in the ",ecoregion, " ecoregion"),
     subtitle =  "1 SE Lambda model")  +
  scale_fill_gradient2(low = "wheat",
                       mid = "darkgreen",
                       high = "red" , 
                       midpoint = 100,   na.value = "lightgrey",
                       limits = c(0,100)) + 
  xlim(st_bbox(plotObs_2)[c(1,3)]) + 
  ylim(st_bbox(plotObs_2)[c(2,4)])

Residuals across the entire temporal extent of the dataset

# rasterize data
plotResid_rast <- pred_glm1 %>% 
         drop_na(resid) %>% 
  #slice_sample(n = 5e4) %>%
  terra::vect(geom = c("Long", "Lat")) %>% 
  terra::set.crs(crs(test_rast)) %>% 
  terra::rasterize(y = test_rast, 
                   field = "resid", 
                   fun = mean) #%>% 
  #terra::aggregate(fact = 2, fun = mean, na.rm = TRUE) %>% 
  #terra::crop(ext(-1950000, 1000000, -1800000, 1000000))

# get the extent of this particular raster, and crop it accordingly
tempExt <- crds(plotResid_rast, na.rm = TRUE)

plotResid_rast_2 <- plotResid_rast %>% 
  crop(ext(min(tempExt[,1]), max(tempExt[,1]),
           min(tempExt[,2]), max(tempExt[,2])) 
       )

# identify locations where residuals are >100 or < -100
badResids_high <- pred_glm1 %>% 
  filter(resid > 100)  %>% 
  terra::vect(geom = c("Long", "Lat")) %>% 
  terra::set.crs(crs(test_rast)) 
badResids_low <- pred_glm1 %>% 
  filter(resid < -100)  %>% 
  terra::vect(geom = c("Long", "Lat")) %>% 
  terra::set.crs(crs(test_rast)) 
# make figures
map <- ggplot() +
geom_spatraster(data =plotResid_rast_2) + 
  geom_sf(data = mapRegions, fill = NA, col = "rosybrown4", lwd = .5) +
  geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_2)),fill=NA )  + 
  geom_sf(data = badResids_high, col = "blue") +
  geom_sf(data = badResids_low, col = "red") +
labs(title = paste0("Resids. (obs. - pred.) from Grass/shrub ecoregion-wide model of ", response),
     subtitle = "bestLambda model \n red points indicate locations that have residuals below -100 \n blue points indicate locatiosn that have residuals above 100") +
  scale_fill_gradient2(low = "red",
                       mid = "white" ,
                       high = "blue" , 
                       midpoint = 0,   na.value = "lightgrey",
                       limits = c(-100,100)
                       ) + 
  xlim(st_bbox(plotObs_2)[c(1,3)]) + 
  ylim(st_bbox(plotObs_2)[c(2,4)])
hist <- ggplot(pred_glm1) + 
  geom_histogram(aes(resid), fill = "lightgrey", col = "darkgrey") + 
  geom_text(aes(x = min(resid)*.9, y = 1500, label = paste0("min = ", round(min(resid),2)))) +
  geom_text(aes(x = max(resid)*.9, y = 1500, label = paste0("max = ", round(max(resid),2))))

library(ggpubr)
ggarrange(map, hist, heights = c(3,1), ncol = 1)

# rasterize data
plotResid_rast <- pred_glm1_1se %>% 
         drop_na(resid) %>% 
  #slice_sample(n = 5e4) %>%
  terra::vect(geom = c("Long", "Lat")) %>% 
  terra::set.crs(crs(test_rast)) %>% 
  terra::rasterize(y = test_rast, 
                   field = "resid", 
                   fun = mean) #%>% 
  #terra::aggregate(fact = 2, fun = mean, na.rm = TRUE) %>% 
  #terra::crop(ext(-1950000, 1000000, -1800000, 1000000))

# get the extent of this particular raster, and crop it accordingly
tempExt <- crds(plotResid_rast, na.rm = TRUE)

plotResid_rast_2 <- plotResid_rast %>% 
  crop(ext(min(tempExt[,1]), max(tempExt[,1]),
           min(tempExt[,2]), max(tempExt[,2])) 
       )

# identify locations where residuals are >100 or < -100
badResids_high <- pred_glm1_1se %>% 
  filter(resid > 100)  %>% 
  terra::vect(geom = c("Long", "Lat")) %>% 
  terra::set.crs(crs(test_rast)) 
badResids_low <- pred_glm1_1se %>% 
  filter(resid < -100)  %>% 
  terra::vect(geom = c("Long", "Lat")) %>% 
  terra::set.crs(crs(test_rast)) 
# make figures
map <- ggplot() +
geom_spatraster(data =plotResid_rast_2) + 
  geom_sf(data = mapRegions, fill = NA, col = "rosybrown4", lwd = .5) +
  geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_2)),fill=NA )  + 
  geom_sf(data = badResids_high, col = "blue") +
  geom_sf(data = badResids_low, col = "red") +
labs(title = paste0("Resids. (obs. - pred.) from Grass/shrub ecoregion-wide model of ", response),
     subtitle = "1 SE Lambda model \n red points indicate locations that have residuals below -100 \n blue points indicate locatiosn that have residuals above 100") +
  scale_fill_gradient2(low = "red",
                       mid = "white" ,
                       high = "blue" , 
                       midpoint = 0,   na.value = "lightgrey",
                       limits = c(-100,100)
                       ) + 
  xlim(st_bbox(plotObs_2)[c(1,3)]) + 
  ylim(st_bbox(plotObs_2)[c(2,4)])
hist <- ggplot(pred_glm1_1se) + 
  geom_histogram(aes(resid), fill = "lightgrey", col = "darkgrey") + 
  geom_text(aes(x = min(resid)*.9, y = 1500, label = paste0("min = ", round(min(resid),2)))) +
  geom_text(aes(x = max(resid)*.9, y = 1500, label = paste0("max = ", round(max(resid),2))))

ggarrange(map, hist, heights = c(3,1), ncol = 1)

### Are there biases of the model predictions across year/lat/long?

# plot residuals against Year
yearResidMod_bestLambda <- ggplot(pred_glm1) + 
  geom_point(aes(x = jitter(Year), y = resid), alpha = .1) + 
  geom_smooth(aes(x = Year, y = resid)) + 
  xlab("Year") + 
  ylab("Residual from best lambda model") +
  ggtitle("from best lamba model")
yearResidMod_1seLambda <- ggplot(pred_glm1_1se) + 
  geom_point(aes(x = jitter(Year), y = resid), alpha = .1) + 
  geom_smooth(aes(x = Year, y = resid)) + 
  xlab("Year") + 
  ylab("Residual from 1 SE lambda model")+
  ggtitle("from 1 SE lamba model")

# plot residuals against Lat
latResidMod_bestLambda <- ggplot(pred_glm1) + 
  geom_point(aes(x = Lat, y = resid), alpha = .1) + 
  geom_smooth(aes(x = Lat, y = resid)) + 
  xlab("Latitude") + 
  ylab("Residual from best lambda model") +
  ggtitle("from best lamba model")
latResidMod_1seLambda <- ggplot(pred_glm1_1se) + 
  geom_point(aes(x = Lat, y = resid), alpha = .1) + 
  geom_smooth(aes(x = Lat, y = resid)) + 
  xlab("Latitude") + 
  ylab("Residual from 1 SE lambda model") +
  ggtitle("from 1 SE lamba model")

# plot residuals against Long
longResidMod_bestLambda <- ggplot(pred_glm1) + 
  geom_point(aes(x = Long, y = resid), alpha = .1) + 
  geom_smooth(aes(x = Long, y = resid)) + 
  xlab("Longitude") + 
  ylab("Residual from best lambda model") +
  ggtitle("from best lamba model")
longResidMod_1seLambda <- ggplot(pred_glm1_1se) + 
  geom_point(aes(x = Long, y = resid), alpha = .1) + 
  geom_smooth(aes(x = Long, y = resid)) + 
  xlab("Longitude") + 
  ylab("Residual from 1 SE lambda model") +
  ggtitle("from 1 SE lamba model")

library(patchwork)
(yearResidMod_bestLambda + yearResidMod_1seLambda) / 
(  latResidMod_bestLambda + latResidMod_1seLambda) /
(  longResidMod_bestLambda + longResidMod_1seLambda)

Quantile plots

Binning predictor variables into “Deciles” (actually percentiles) and looking at the mean predicted probability for each percentile. The use of the word Decentiles is just a legacy thing (they started out being actual Percentiles)

# get deciles for best lambda model 
if (length(prednames_fig) == 0) {
  print("The best lambda model only contains one predictor (an intercept), so decile plots aren't possible to generate")
} else {
  pred_glm1_deciles <- predvars2deciles(pred_glm1,
                                      response_vars = response_vars,
                                        pred_vars = prednames_fig, 
                                       cut_points = seq(0, 1, 0.005))
}
# get deciles for 1 SE lambda model 
if (length(prednames_fig_1se) == 0) {
  print("The 1SE lambda model only contains one predictor (an intercept), so decile plots aren't possible to generate")
} else {
  pred_glm1_deciles_1se <- predvars2deciles(pred_glm1_1se,
                                      response_vars = response_vars,
                                        pred_vars = prednames_fig_1se, 
                                       cut_points = seq(0, 1, 0.005))
}
## [1] "The 1SE lambda model only contains one predictor (an intercept), so decile plots aren't possible to generate"

Here is a quick version of LOESS curves fit to raw data (to double-check the quantile plot calculations)

if (length(prednames_fig) == 0) {
  print("The model only contains one predictor (an intercept), so decile plots aren't possible to generate")
} else {
  pred_glm1 %>%
  select(all_of(c(prednames_fig, response_vars))) %>%
  pivot_longer(cols = prednames_fig)  %>%
  ggplot() +
  facet_wrap(~name, scales = "free") +
  geom_point(aes(x = value, y =  .data[[paste(response)]]), col = "darkblue", alpha = .1)  + # observed values
  geom_point(aes(x = value, y = .data[[response_vars[2]]]), col = "lightblue", alpha = .1) + # model-predicted values
  geom_smooth(aes(x = value, y =  .data[[paste(response)]]), col = "black", se = FALSE) +
  geom_smooth(aes(x = value, y = .data[[response_vars[2]]]), col = "brown", se = FALSE)

}

Below are the actual quantile plots (note that the predictor variables are scaled)

if (length(prednames_fig) == 0) {
  print("The model only contains one predictor (an intercept), so decile plots aren't possible to generate")
} else {

# publication quality version
g3 <- decile_dotplot_pq(pred_glm1_deciles, response = response, IQR = TRUE) + ggtitle("Decile Plot")

g4 <- add_dotplot_inset(g3, df = pred_glm1_deciles, dfRaw = pred_glm1, add_smooth = TRUE, deciles = FALSE)

  
if(save_figs) {
  png(paste0("figures/quantile_plots/quantile_plot_", response,  "_",ecoregion,".png"), 
     units = "in", res = 600, width = 5.5, height = 3.5 )
    print(g4)
  dev.off()
}

g4
}

if (length(prednames_fig_1se) == 0) {
  print("The 1 se lambda model only contains one predictor (an intercept), so decile plots aren't possible to generate")

  } else {

# publication quality version
g3 <- decile_dotplot_pq(pred_glm1_deciles_1se, response = response, IQR = TRUE) + ggtitle("Decile Plot")

g4 <- add_dotplot_inset(g3, df = pred_glm1_deciles_1se, dfRaw = pred_glm1_1se, add_smooth = TRUE, deciles = FALSE)

  
if(save_figs) {
  png(paste0("figures/quantile_plots/quantile_plot_", response,  "_",ecoregion,".png"), 
     units = "in", res = 600, width = 5.5, height = 3.5 )
    print(g4)
  dev.off()
}

g4
}
## [1] "The 1 se lambda model only contains one predictor (an intercept), so decile plots aren't possible to generate"

Deciles Filtered

20th and 80th percentiles for each climate variable

df <- pred_glm1[, prednames_fig] #%>% 
  #mutate(MAT = MAT - 273.15) # k to c
quantiles <- map(df, quantile, probs = c(0.2, 0.8), na.rm = TRUE)

Filtered ‘Decile’ plots of data. These plots show each vegetation variable, but only based on data that falls into the upper and lower two deciles of each predictor variable.

if (length(prednames_fig) == 0) {
  print("The model only contains one predictor (an intercept), so decile plots aren't possible to generate")
} else {
pred_glm1_deciles_filt <- predvars2deciles( pred_glm1, 
                         response_vars = response_vars,
                         pred_vars = prednames_fig,
                         filter_var = TRUE,
                         filter_vars = prednames_fig,
                         cut_points = seq(0, 1, 0.005)) 

decile_dotplot_filtered_pq(pred_glm1_deciles_filt, xvars = prednames_fig)
#decile_dotplot_filtered_pq(pred_glm1_deciles_filt)

}
## Processed 16460 groups out of 217362. 8% done. Time elapsed: 3s. ETA: 36s.Processed 22103 groups out of 217362. 10% done. Time elapsed: 4s. ETA: 35s.Processed 27764 groups out of 217362. 13% done. Time elapsed: 5s. ETA: 34s.Processed 32355 groups out of 217362. 15% done. Time elapsed: 6s. ETA: 34s.Processed 37747 groups out of 217362. 17% done. Time elapsed: 7s. ETA: 33s.Processed 43206 groups out of 217362. 20% done. Time elapsed: 8s. ETA: 32s.Processed 48783 groups out of 217362. 22% done. Time elapsed: 9s. ETA: 31s.Processed 54330 groups out of 217362. 25% done. Time elapsed: 10s. ETA: 30s.Processed 59788 groups out of 217362. 28% done. Time elapsed: 11s. ETA: 29s.Processed 65346 groups out of 217362. 30% done. Time elapsed: 12s. ETA: 28s.Processed 70920 groups out of 217362. 33% done. Time elapsed: 13s. ETA: 26s.Processed 76471 groups out of 217362. 35% done. Time elapsed: 14s. ETA: 25s.Processed 82053 groups out of 217362. 38% done. Time elapsed: 15s. ETA: 24s.Processed 87027 groups out of 217362. 40% done. Time elapsed: 16s. ETA: 24s.Processed 92625 groups out of 217362. 43% done. Time elapsed: 17s. ETA: 22s.Processed 98286 groups out of 217362. 45% done. Time elapsed: 18s. ETA: 21s.Processed 103935 groups out of 217362. 48% done. Time elapsed: 19s. ETA: 20s.Processed 109439 groups out of 217362. 50% done. Time elapsed: 20s. ETA: 19s.Processed 115010 groups out of 217362. 53% done. Time elapsed: 21s. ETA: 18s.Processed 120548 groups out of 217362. 55% done. Time elapsed: 22s. ETA: 17s.Processed 126040 groups out of 217362. 58% done. Time elapsed: 23s. ETA: 16s.Processed 131580 groups out of 217362. 61% done. Time elapsed: 24s. ETA: 15s.Processed 137082 groups out of 217362. 63% done. Time elapsed: 25s. ETA: 14s.Processed 142575 groups out of 217362. 66% done. Time elapsed: 26s. ETA: 13s.Processed 147912 groups out of 217362. 68% done. Time elapsed: 27s. ETA: 12s.Processed 153477 groups out of 217362. 71% done. Time elapsed: 28s. ETA: 11s.Processed 159108 groups out of 217362. 73% done. Time elapsed: 29s. ETA: 10s.Processed 164561 groups out of 217362. 76% done. Time elapsed: 30s. ETA: 9s.Processed 169742 groups out of 217362. 78% done. Time elapsed: 31s. ETA: 8s.Processed 175251 groups out of 217362. 81% done. Time elapsed: 32s. ETA: 7s.Processed 180761 groups out of 217362. 83% done. Time elapsed: 33s. ETA: 6s.Processed 186341 groups out of 217362. 86% done. Time elapsed: 34s. ETA: 5s.Processed 191938 groups out of 217362. 88% done. Time elapsed: 35s. ETA: 4s.Processed 197562 groups out of 217362. 91% done. Time elapsed: 36s. ETA: 3s.Processed 203183 groups out of 217362. 93% done. Time elapsed: 37s. ETA: 2s.Processed 208721 groups out of 217362. 96% done. Time elapsed: 38s. ETA: 1s.Processed 214340 groups out of 217362. 99% done. Time elapsed: 39s. ETA: 0s.Processed 217362 groups out of 217362. 100% done. Time elapsed: 39s. ETA: 0s.

Filtered quantile figure with middle 2 deciles also shown

if (length(prednames_fig) == 0) {
  print("The model only contains one predictor (an intercept), so decile plots aren't possible to generate")
} else {
pred_glm1_deciles_filt_mid <- predvars2deciles(pred_glm1, 
                         response_vars = response_vars,
                         pred_vars = prednames_fig,
                         filter_vars = prednames_fig,
                         filter_var = TRUE,
                         add_mid = TRUE,
                         cut_points = seq(0, 1, 0.005))

g <- decile_dotplot_filtered_pq(df = pred_glm1_deciles_filt_mid, xvars = prednames_fig)
g

if(save_figs) {x
jpeg(paste0("figures/quantile_plots/quantile_plot_filtered_mid_v1", , ".jpeg"),
     units = "in", res = 600, width = 5.5, height = 6 )
  g 
dev.off()
}
}
## Processed 15667 groups out of 326145. 5% done. Time elapsed: 3s. ETA: 59s.Processed 21060 groups out of 326145. 6% done. Time elapsed: 4s. ETA: 57s.Processed 26437 groups out of 326145. 8% done. Time elapsed: 5s. ETA: 56s.Processed 31768 groups out of 326145. 10% done. Time elapsed: 6s. ETA: 55s.Processed 37031 groups out of 326145. 11% done. Time elapsed: 7s. ETA: 54s.Processed 41341 groups out of 326145. 13% done. Time elapsed: 8s. ETA: 55s.Processed 46310 groups out of 326145. 14% done. Time elapsed: 9s. ETA: 54s.Processed 51345 groups out of 326145. 16% done. Time elapsed: 10s. ETA: 53s.Processed 56521 groups out of 326145. 17% done. Time elapsed: 11s. ETA: 52s.Processed 61951 groups out of 326145. 19% done. Time elapsed: 12s. ETA: 51s.Processed 67333 groups out of 326145. 21% done. Time elapsed: 13s. ETA: 49s.Processed 72774 groups out of 326145. 22% done. Time elapsed: 14s. ETA: 48s.Processed 78210 groups out of 326145. 24% done. Time elapsed: 15s. ETA: 47s.Processed 83628 groups out of 326145. 26% done. Time elapsed: 16s. ETA: 46s.Processed 88947 groups out of 326145. 27% done. Time elapsed: 17s. ETA: 45s.Processed 94329 groups out of 326145. 29% done. Time elapsed: 18s. ETA: 44s.Processed 99673 groups out of 326145. 31% done. Time elapsed: 19s. ETA: 43s.Processed 105087 groups out of 326145. 32% done. Time elapsed: 20s. ETA: 42s.Processed 110497 groups out of 326145. 34% done. Time elapsed: 21s. ETA: 41s.Processed 115929 groups out of 326145. 36% done. Time elapsed: 22s. ETA: 40s.Processed 121071 groups out of 326145. 37% done. Time elapsed: 23s. ETA: 39s.Processed 126210 groups out of 326145. 39% done. Time elapsed: 24s. ETA: 38s.Processed 131149 groups out of 326145. 40% done. Time elapsed: 25s. ETA: 37s.Processed 136426 groups out of 326145. 42% done. Time elapsed: 26s. ETA: 36s.Processed 141704 groups out of 326145. 43% done. Time elapsed: 27s. ETA: 35s.Processed 146880 groups out of 326145. 45% done. Time elapsed: 28s. ETA: 34s.Processed 152289 groups out of 326145. 47% done. Time elapsed: 29s. ETA: 33s.Processed 157691 groups out of 326145. 48% done. Time elapsed: 30s. ETA: 32s.Processed 163071 groups out of 326145. 50% done. Time elapsed: 31s. ETA: 31s.Processed 168438 groups out of 326145. 52% done. Time elapsed: 32s. ETA: 30s.Processed 173785 groups out of 326145. 53% done. Time elapsed: 33s. ETA: 29s.Processed 179133 groups out of 326145. 55% done. Time elapsed: 34s. ETA: 28s.Processed 184543 groups out of 326145. 57% done. Time elapsed: 35s. ETA: 27s.Processed 189959 groups out of 326145. 58% done. Time elapsed: 36s. ETA: 25s.Processed 195340 groups out of 326145. 60% done. Time elapsed: 37s. ETA: 24s.Processed 200725 groups out of 326145. 62% done. Time elapsed: 38s. ETA: 23s.Processed 206131 groups out of 326145. 63% done. Time elapsed: 39s. ETA: 22s.Processed 211458 groups out of 326145. 65% done. Time elapsed: 40s. ETA: 21s.Processed 216856 groups out of 326145. 66% done. Time elapsed: 41s. ETA: 20s.Processed 222064 groups out of 326145. 68% done. Time elapsed: 42s. ETA: 19s.Processed 227487 groups out of 326145. 70% done. Time elapsed: 43s. ETA: 18s.Processed 232879 groups out of 326145. 71% done. Time elapsed: 44s. ETA: 17s.Processed 238274 groups out of 326145. 73% done. Time elapsed: 45s. ETA: 16s.Processed 243680 groups out of 326145. 75% done. Time elapsed: 46s. ETA: 15s.Processed 248623 groups out of 326145. 76% done. Time elapsed: 47s. ETA: 14s.Processed 249846 groups out of 326145. 77% done. Time elapsed: 48s. ETA: 14s.Processed 255023 groups out of 326145. 78% done. Time elapsed: 49s. ETA: 13s.Processed 260619 groups out of 326145. 80% done. Time elapsed: 50s. ETA: 12s.Processed 266230 groups out of 326145. 82% done. Time elapsed: 51s. ETA: 11s.Processed 271798 groups out of 326145. 83% done. Time elapsed: 52s. ETA: 10s.Processed 277453 groups out of 326145. 85% done. Time elapsed: 53s. ETA: 9s.Processed 282709 groups out of 326145. 87% done. Time elapsed: 54s. ETA: 8s.Processed 288337 groups out of 326145. 88% done. Time elapsed: 55s. ETA: 7s.Processed 293988 groups out of 326145. 90% done. Time elapsed: 56s. ETA: 6s.Processed 299616 groups out of 326145. 92% done. Time elapsed: 57s. ETA: 5s.Processed 305076 groups out of 326145. 94% done. Time elapsed: 58s. ETA: 4s.Processed 310547 groups out of 326145. 95% done. Time elapsed: 59s. ETA: 2s.Processed 316202 groups out of 326145. 97% done. Time elapsed: 60s. ETA: 1s.Processed 321808 groups out of 326145. 99% done. Time elapsed: 61s. ETA: 0s.Processed 326145 groups out of 326145. 100% done. Time elapsed: 62s. ETA: 0s.

Cross-validation

Using best lambda model

Use terms from global model to re-fit and predict on different held out regions

Figures show residuals for each of the models fit to held-out ecoregions

These models were fit to six ecoregions, and then predict on the indicated heldout ecoregion

if (length(prednames_fig) == 0) {
  print("The model only contains one predictor (an intercept), so cross validation isn't practical")
} else {
  
## code from Tredennick et al. 2020
# try each separate level II ecoregion as a test set
# make a list to hold output data
outList <- vector(mode = "list", length = length(sort(unique(modDat_1$NA_L2NAME))))
# obs_pred <- data.frame(ecoregion = character(),obs = numeric(),
#                        pred_opt = numeric(), pred_null = numeric()#,
#                        #pred_nopenalty = numeric()
#                        )

## get the model specification from the global model
mat <- as.matrix(coef(fit_glm_bestLambda, s = "lambda.min"))
mat2 <- mat[mat[,1] != 0,]

f_cv <- as.formula(paste0(response, " ~ ", paste0(names(mat2)[2:length(names(mat2))], collapse = " + ")))

X_cv <- model.matrix(object = f_cv, data = modDat_1_s)
# get response variable
y_cv <- as.matrix(modDat_1_s[,response])

  
# now, loop through so with each iteration, a different ecoregion is held out
 for(i_eco in sort(unique(modDat_1_s$NA_L2NAME))){

  # split into training and test sets
  test_eco <- i_eco
  print(test_eco)
  # identify the rowID of observations to be in the training and test datasets
  train <- which(modDat_1_s$NA_L2NAME!=test_eco) # data for all ecoregions that aren't 'i_eco'
  test <- which(modDat_1_s$NA_L2NAME==test_eco) # data for the ecoregion that is 'i_eco'

  trainDat_all <- modDat_1_s %>% 
    slice(train) %>% 
    select(-newRegion)
  testDat_all <- modDat_1_s %>% 
    slice(test) %>% 
    select(-newRegion)

  # get the model matrices for input and response variables for cross validation model specification
  X_train <- as.matrix(X_cv[train,])
  X_test <- as.matrix(X_cv[test,])

  y_train <- modDat_1_s[train,response]
  y_test <- modDat_1_s[test,response]
  
  # get the model matrices for input and response variables for original model specification
  X_train_glob <- as.matrix(X[train,])
  X_test_glob <- as.matrix(X[test,])

  y_train_glob <- modDat_1_s[train,response]
  y_test_glob <- modDat_1_s[test,response]

  train_eco <- modDat_1_s$NA_L2NAME[train]

  ## just try a regular glm w/ the components from the global model
  fit_i <- glm(data = trainDat_all, formula = f_cv, 
    ,
               family =  stats::Gamma(link = "log")
    )
    
  # lasso model predictions with the optimal lambda
  optimal_pred <- predict(fit_i, newdata= testDat_all, type = "response"
                          )
  # null model and predictions
  # the "null" model in this case is the global model
  # predict on the test data for this iteration w/ the global model 
  null_pred <- predict.glm(fit_glm_bestLambda, newdata = testDat_all, type = "response")

  
  # save data
  tmp <- data.frame(ecoRegion_holdout = rep(test_eco,length(y_test)),obs=y_test,
                    pred_opt=optimal_pred, pred_null=null_pred#,
                    #pred_nopenalty=nopen_pred
                    ) %>%
    cbind(testDat_all)
  
  # calculate RMSE, bias, etc. of 
  # RMSE of CV model 
  RMSE_optimal <- yardstick::rmse(data = data.frame(optimal_pred, y_test), truth = "y_test", estimate = "optimal_pred")[1,]$.estimate
  # RMSE of global model
  RMSE_null <- yardstick::rmse(data = data.frame(null_pred, y_test), truth = "y_test", estimate = "null_pred")[1,]$.estimate
  # bias of CV model
  bias_optimal <- mean(y_test - optimal_pred)
  # bias of global model
  bias_null <-  mean( y_test - null_pred )
  
  # put output into a list
  tmpList <- list("testRegion" = i_eco,
    "modelObject" = fit_i,
       "modelPredictions" = tmp, 
    "performanceMetrics" = data.frame("RMSE_cvModel" = RMSE_optimal, 
                                      "RMSE_globalModel" = RMSE_null, 
                                      "bias_cvModel" = bias_optimal, 
                                      "bias_globalModel" = bias_null))

  # save model outputs
  outList[[which(sort(unique(modDat_1_s$NA_L2NAME)) == i_eco)]] <- tmpList
 }
}
## [1] "COLD DESERTS"
## [1] "MEDITERRANEAN PIEDMONT"
## [1] "SEMIARID PLAIN AND PRAIRIES"
## [1] "TEMPERATE PRAIRIES"
## [1] "WARM DESERTS"
## [1] "WEST-CENTRAL SEMIARID PRAIRIES"

Below are the RMSE and bias values for predictions made for each holdout level II ecoregion, compared to predictions from the global model for that same ecoregion

# table of model performance
map(outList, .f = function(x) {
  cbind(data.frame("holdout region" = x$testRegion),  x$performanceMetrics)
}
) %>% 
  purrr::list_rbind() %>% 
  kable(col.names = c("Held-out ecoregion", "RMSE of CV model", "RMSE of global model", 
                      "bias of CV model - mean(obs-pred.)", "bias of global model- mean(obs-pred.)"), 
        caption = "Performance of Cross Validation using 'best lambda' model specification") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed")) 
Performance of Cross Validation using ‘best lambda’ model specification
Held-out ecoregion RMSE of CV model RMSE of global model bias of CV model - mean(obs-pred.) bias of global model- mean(obs-pred.)
COLD DESERTS 26.88670 17.42212 -17.398240 -0.7546112
MEDITERRANEAN PIEDMONT 23.49523 20.51927 -3.935107 -0.3336532
SEMIARID PLAIN AND PRAIRIES 29.79030 24.68797 7.285292 -0.2116874
TEMPERATE PRAIRIES 64.64083 39.45014 -37.271386 -14.6141885
WARM DESERTS 16.41719 15.79510 -1.893220 -0.7810884
WEST-CENTRAL SEMIARID PRAIRIES 26.07172 23.03745 13.088519 4.6672097
# visualize model predictions
for (i in 1:length(unique(modDat_1_s$NA_L2NAME))) {
  holdoutRegion <- outList[[i]]$testRegion
  predictionData <- outList[[i]]$modelPredictions
  modTerms <- as.matrix(coef(outList[[i]]$modelObject)) %>%
    as.data.frame() %>%
    filter(V1!=0) %>%
    rownames()

  # calculate residuals
  predictionData <- predictionData %>%
  mutate(resid = .[["obs"]] - .[["pred_opt"]] ,
         resid_globMod = .[["obs"]]  - .[["pred_null"]])


# rasterize
# use 'test_rast' from earlier

  # rasterize data
plotObs <- predictionData %>%
         drop_na(paste(response)) %>%
  #slice_sample(n = 5e4) %>%
  terra::vect(geom = c("Long", "Lat")) %>%
  terra::set.crs(crs(test_rast)) %>%
  terra::rasterize(y = test_rast,
                   field = "resid",
                   fun = mean) #%>%
  #terra::aggregate(fact = 2, fun = mean, na.rm = TRUE) %>%
  #terra::crop(ext(-1950000, 1000000, -1800000, 1000000))

tempExt <- crds(plotObs, na.rm = TRUE)

plotObs_2 <- plotObs %>% 
  crop(ext(min(tempExt[,1]), max(tempExt[,1]),
           min(tempExt[,2]), max(tempExt[,2])) 
       )

# identify locations where residuals are >100 or < -100
badResids_high <- predictionData %>% 
  filter(resid > 100)  %>% 
  terra::vect(geom = c("Long", "Lat")) %>% 
  terra::set.crs(crs(test_rast)) 
badResids_low <- predictionData %>% 
  filter(resid < -100)  %>% 
  terra::vect(geom = c("Long", "Lat")) %>% 
  terra::set.crs(crs(test_rast)) 


# make figures
# make histogram
hist_i <- ggplot(predictionData) +
  geom_histogram(aes(resid_globMod), col = "darkgrey", fill = "lightgrey") +
  xlab(c("Residuals (obs. - pred.)"))
# make map
map_i <-  ggplot() +
geom_spatraster(data = plotObs_2) +
  geom_sf(data = mapRegions, fill = NA, col = "rosybrown4", lwd = .5) +
  geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_2)),fill=NA ) +
  geom_sf(data = badResids_high, col = "blue") +
  geom_sf(data = badResids_low, col = "red") +
labs(title = paste0("Residuals (obs. - pred.) for predictions of \n", holdoutRegion, " \n from a model fit to other ecoregions"),
     subtitle = paste0(response, " ~ ", paste0( modTerms, collapse = " + "))) +
  scale_fill_gradient2(low = "red",
                       mid = "white" ,
                       high = "blue" ,
                       midpoint = 0,   na.value = "lightgrey",
                       limits = c(-100, 100))  + 
  xlim(st_bbox(plotObs_2)[c(1,3)]) + 
  ylim(st_bbox(plotObs_2)[c(2,4)])

 assign(paste0("residPlot_",holdoutRegion),
   value = ggarrange(map_i, hist_i, heights = c(3,1), ncol = 1)
)

}

  lapply(unique(modDat_1_s$NA_L2NAME), FUN = function(x) {
    get(paste0("residPlot_", x))
  })
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

Using 1se lambda model

Use terms from global model to re-fit and predict on different held out regions

Figures show residuals for each of the models fit to held-out ecoregions

These models were fit to six ecoregions, and then predict on the indicated heldout ecoregion

if (length(prednames_fig_1se) == 0) {
  print("The model only contains one predictor (an intercept), so cross validation isn't practical")
} else {

## code from Tredennick et al. 2020
# try each separate level II ecoregion as a test set
# make a list to hold output data
outList <- vector(mode = "list", length = length(sort(unique(modDat_1$NA_L2NAME))))
# obs_pred <- data.frame(ecoregion = character(),obs = numeric(),
#                        pred_opt = numeric(), pred_null = numeric()#,
#                        #pred_nopenalty = numeric()
#                        )

## get the model specification from the global model
mat <- as.matrix(coef(fit_glm_se, s = "lambda.min"))
mat2 <- mat[mat[,1] != 0,]

f_cv <- as.formula(paste0(response, " ~ ", paste0(names(mat2)[2:length(names(mat2))], collapse = " + ")))

X_cv <- model.matrix(object = f_cv, data = modDat_1_s)
# get response variable
y_cv <- as.matrix(modDat_1_s[,response])

  
# now, loop through so with each iteration, a different ecoregion is held out
 for(i_eco in sort(unique(modDat_1_s$NA_L2NAME))){

  # split into training and test sets
  test_eco <- i_eco
  print(test_eco)
  # identify the rowID of observations to be in the training and test datasets
  train <- which(modDat_1_s$NA_L2NAME!=test_eco) # data for all ecoregions that aren't 'i_eco'
  test <- which(modDat_1_s$NA_L2NAME==test_eco) # data for the ecoregion that is 'i_eco'

  trainDat_all <- modDat_1_s %>% 
    slice(train) %>% 
    select(-newRegion)
  testDat_all <- modDat_1_s %>% 
    slice(test) %>% 
    select(-newRegion)

  # get the model matrices for input and response variables for cross validation model specification
  X_train <- as.matrix(X_cv[train,])
  X_test <- as.matrix(X_cv[test,])

  y_train <- modDat_1_s[train,response]
  y_test <- modDat_1_s[test,response]
  
  # get the model matrices for input and response variables for original model specification
  X_train_glob <- as.matrix(X[train,])
  X_test_glob <- as.matrix(X[test,])

  y_train_glob <- modDat_1_s[train,response]
  y_test_glob <- modDat_1_s[test,response]

  train_eco <- modDat_1_s$NA_L2NAME[train]

  ## just try a regular glm w/ the components from the global model
  fit_i <- glm(data = trainDat_all, formula = f_cv, 
    ,
               family =  stats::Gamma(link = "log")
    )

    coef(fit_i)
    
  # lasso model predictions with the optimal lambda
  optimal_pred <- predict(fit_i, newdata= testDat_all, type = "response"
                          )
  # null model and predictions
  # the "null" model in this case is the global model
  # predict on the test data for this iteration w/ the global model 
  null_pred <- predict.glm(fit_glm_se, newdata = testDat_all, type = "response")

  # save data
  tmp <- data.frame(ecoRegion_holdout = rep(test_eco,length(y_test)),obs=y_test,
                    pred_opt=optimal_pred, pred_null=null_pred#,
                    #pred_nopenalty=nopen_pred
                    ) %>%
    cbind(testDat_all)
    
  # calculate RMSE, bias, etc. of 
  # RMSE of CV model 
  RMSE_optimal <- yardstick::rmse(data = data.frame(optimal_pred, y_test), truth = "y_test", estimate = "optimal_pred")[1,]$.estimate
  # RMSE of global model
  RMSE_null <- yardstick::rmse(data = data.frame(null_pred, y_test), truth = "y_test", estimate = "null_pred")[1,]$.estimate
  # bias of CV model
  bias_optimal <- mean(y_test - optimal_pred)
  # bias of global model
  bias_null <-  mean(y_test - null_pred )
  
  # put output into a list
  tmpList <- list("testRegion" = i_eco,
    "modelObject" = fit_i,
       "modelPredictions" = tmp, 
    "performanceMetrics" = data.frame("RMSE_cvModel" = RMSE_optimal, 
                                      "RMSE_globalModel" = RMSE_null, 
                                      "bias_cvModel" = bias_optimal, 
                                      "bias_globalModel" = bias_null))

  # save model outputs
  outList[[which(sort(unique(modDat_1_s$NA_L2NAME)) == i_eco)]] <- tmpList
 }
}
## [1] "The model only contains one predictor (an intercept), so cross validation isn't practical"

Below are the RMSE and bias values for predictions made for each holdout level II ecoregion, compared to predictions from the global model for that same ecoregion

if (length(prednames_fig_1se) == 0) {
  print("The model only contains one predictor (an intercept), so cross validation isn't practical")
} else {
# table of model performance
map(outList, .f = function(x) {
  cbind(data.frame("holdout region" = x$testRegion),  x$performanceMetrics)
}
) %>% 
  purrr::list_rbind() %>% 
  kable(col.names = c("Held-out ecoregion", "RMSE of CV model", "RMSE of global model", 
                      "bias of CV model - mean(obs-pred.)", "bias of global model - mean(obs-pred.)"), 
        caption = "Performance of Cross Validation using '1 SE lambda' model specification") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed")) 
}
## [1] "The model only contains one predictor (an intercept), so cross validation isn't practical"
if (length(prednames_fig_1se) == 0) {
  print("The model only contains one predictor (an intercept), so cross validation isn't practical")
} else {
for (i in 1:length(unique(modDat_1_s$NA_L2NAME))) {
  holdoutRegion <- outList[[i]]$testRegion
  predictionData <- outList[[i]]$modelPredictions
  modTerms <- as.matrix(coef(outList[[i]]$modelObject)) %>%
    as.data.frame() %>%
    filter(V1!=0) %>%
    rownames()

  # calculate residuals
  predictionData <- predictionData %>%
  mutate(resid = .[["obs"]] - .[["pred_opt"]] ,
         resid_globMod = .[["obs"]]  - .[["pred_null"]])


# rasterize
# use 'test_rast' from earlier

  # rasterize data
plotObs <- predictionData %>%
         drop_na(paste(response)) %>%
  #slice_sample(n = 5e4) %>%
  terra::vect(geom = c("Long", "Lat")) %>%
  terra::set.crs(crs(test_rast)) %>%
  terra::rasterize(y = test_rast,
                   field = "resid",
                   fun = mean) #%>%
  #terra::aggregate(fact = 2, fun = mean, na.rm = TRUE) %>%
  #terra::crop(ext(-1950000, 1000000, -1800000, 1000000))

tempExt <- crds(plotObs, na.rm = TRUE)

plotObs_2 <- plotObs %>% 
  crop(ext(min(tempExt[,1]), max(tempExt[,1]),
           min(tempExt[,2]), max(tempExt[,2])) 
       )

# identify locations where residuals are >100 or < -100
badResids_high <- predictionData %>% 
  filter(resid > 100)  %>% 
  terra::vect(geom = c("Long", "Lat")) %>% 
  terra::set.crs(crs(test_rast)) 
badResids_low <- predictionData %>% 
  filter(resid < -100)  %>% 
  terra::vect(geom = c("Long", "Lat")) %>% 
  terra::set.crs(crs(test_rast)) 


# make figures
# make histogram
hist_i <- ggplot(predictionData) +
  geom_histogram(aes(resid_globMod), col = "darkgrey", fill = "lightgrey") +
  xlab(c("Residuals (obs. - pred.)"))
# make map
map_i <-  ggplot() +
geom_spatraster(data = plotObs_2) +
  geom_sf(data = mapRegions, fill = NA, col = "rosybrown4", lwd = .5) +
  geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_2)),fill=NA ) +
  geom_sf(data = badResids_high, col = "blue") +
  geom_sf(data = badResids_low, col = "red") +
labs(title = paste0("Residuals (obs. - pred.) for predictions of \n", holdoutRegion, " \n from a model fit to other ecoregions"),
     subtitle = paste0(response, " ~ ", paste0( modTerms, collapse = " + "))) +
  scale_fill_gradient2(low = "red",
                       mid = "white" ,
                       high = "blue" ,
                       midpoint = 0,   na.value = "lightgrey",
                       limits = c(-100, 100))  + 
  xlim(st_bbox(plotObs_2)[c(1,3)]) + 
  ylim(st_bbox(plotObs_2)[c(2,4)])

 assign(paste0("residPlot_",holdoutRegion),
   value = ggarrange(map_i, hist_i, heights = c(3,1), ncol = 1)
)

}

  lapply(unique(modDat_1_s$NA_L2NAME), FUN = function(x) {
    get(paste0("residPlot_", x))
  })
}
## [1] "The model only contains one predictor (an intercept), so cross validation isn't practical"

Save output

# # glm models
# #mods2save <- butcher::butcher(mod_glmFinal) # removes some model components so the saved object isn't huge
# 
# #mods2save$formula <- best_form
# #mods2save$pred_vars_inter <- pred_vars_inter # so have interactions
# #n <- nrow(df_sample)
# #mods2save$data_rows <- n
# 
# 
# if(!test_run) {
#   saveRDS(mods2save, 
#         paste0("./models/glm_beta_model_CONUSwide_", s, "_n", n, 
#         #sample_group, 
#         ".RDS"))
#   if (byRegion == TRUE) {
#     ## western forests
#      saveRDS(mods2save_WF, 
#         paste0("./models/glm_beta_model_WesternForests_", s, "_n", n, 
#         #sample_group, 
#         ".RDS"))
#     ## eastern forests
#      saveRDS(mods2save_EF, 
#         paste0("./models/glm_beta_model_EasternForests_", s, "_n", n, 
#         #sample_group, 
#         ".RDS"))
#      ## grass/shrub
#      saveRDS(mods2save_G, 
#         paste0("./models/glm_beta_model_GrassShrub_", s, "_n", n, 
#         #sample_group, 
#         ".RDS"))
#   }
# }
## partial dependence plots
#vip::vip(mod_glmFinal, num_features = 15)

#pdp_all_vars(mod_glmFinal, mod_vars = pred_vars, ylab = 'probability',train = df_small)

#caret::varImp(fit)

session info

Hash of current commit (i.e. to ID the version of the code used)

system("git rev-parse HEAD", intern=TRUE)
## [1] "79890c55a196d40eb16ae968701c4515b44c260c"

Packages etc.

sessionInfo()
## R version 4.4.0 (2024-04-24)
## Platform: aarch64-apple-darwin20
## Running under: macOS Sonoma 14.7.5
## 
## Matrix products: default
## BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/Denver
## tzcode source: internal
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] doMC_1.3.8                 iterators_1.0.14           foreach_1.5.2              ggpubr_0.6.0               factoextra_1.0.7          
##  [6] USA.state.boundaries_1.0.1 glmnet_4.1-8               Matrix_1.7-0               kableExtra_1.4.0           rsample_1.2.1             
## [11] here_1.0.1                 StepBeta_2.1.0             ggtext_0.1.2               knitr_1.49                 gridExtra_2.3             
## [16] pdp_0.8.2                  GGally_2.2.1               lubridate_1.9.4            forcats_1.0.0              stringr_1.5.1             
## [21] dplyr_1.1.4                purrr_1.0.4                readr_2.1.5                tidyr_1.3.1                tibble_3.2.1              
## [26] tidyverse_2.0.0            caret_6.0-94               lattice_0.22-6             ggplot2_3.5.1              sf_1.0-20                 
## [31] tidyterra_0.6.1            terra_1.8-21               ggspatial_1.1.9            dtplyr_1.3.1               patchwork_1.3.0           
## 
## loaded via a namespace (and not attached):
##   [1] RColorBrewer_1.1-3   rstudioapi_0.17.1    jsonlite_1.9.1       shape_1.4.6.1        magrittr_2.0.3       modeltools_0.2-23   
##   [7] farver_2.1.2         rmarkdown_2.29       vctrs_0.6.5          rstatix_0.7.2        htmltools_0.5.8.1    broom_1.0.7         
##  [13] Formula_1.2-5        pROC_1.18.5          sass_0.4.9           parallelly_1.37.1    KernSmooth_2.23-22   bslib_0.9.0         
##  [19] plyr_1.8.9           sandwich_3.1-0       zoo_1.8-12           cachem_1.1.0         commonmark_1.9.1     lifecycle_1.0.4     
##  [25] pkgconfig_2.0.3      R6_2.6.1             fastmap_1.2.0        future_1.33.2        digest_0.6.37        colorspace_2.1-1    
##  [31] furrr_0.3.1          rprojroot_2.0.4      pkgload_1.3.4        labeling_0.4.3       yardstick_1.3.1      timechange_0.3.0    
##  [37] mgcv_1.9-1           abind_1.4-8          compiler_4.4.0       proxy_0.4-27         aod_1.3.3            withr_3.0.2         
##  [43] backports_1.5.0      carData_3.0-5        betareg_3.1-4        DBI_1.2.3            ggstats_0.9.0        ggsignif_0.6.4      
##  [49] MASS_7.3-60.2        lava_1.8.0           classInt_0.4-10      gtools_3.9.5         ModelMetrics_1.2.2.2 tools_4.4.0         
##  [55] units_0.8-5          lmtest_0.9-40        future.apply_1.11.2  nnet_7.3-19          glue_1.8.0           nlme_3.1-164        
##  [61] gridtext_0.1.5       grid_4.4.0           reshape2_1.4.4       generics_0.1.3       recipes_1.1.0        gtable_0.3.6        
##  [67] tzdb_0.4.0           class_7.3-22         data.table_1.17.0    hms_1.1.3            utf8_1.2.4           car_3.1-2           
##  [73] xml2_1.3.7           flexmix_2.3-19       markdown_1.13        ggrepel_0.9.5        pillar_1.10.1        splines_4.4.0       
##  [79] survival_3.5-8       tidyselect_1.2.1     svglite_2.1.3        stats4_4.4.0         xfun_0.51            hardhat_1.4.0       
##  [85] timeDate_4032.109    stringi_1.8.4        yaml_2.3.10          evaluate_1.0.3       codetools_0.2-20     cli_3.6.4           
##  [91] rpart_4.1.23         systemfonts_1.2.1    munsell_0.5.1        jquerylib_0.1.4      Rcpp_1.0.14          globals_0.16.3      
##  [97] gower_1.0.1          listenv_0.9.1        viridisLite_0.4.2    ipred_0.9-15         scales_1.3.0         prodlim_2024.06.25  
## [103] e1071_1.7-14         crayon_1.5.3         combinat_0.0-8       rlang_1.1.5          cowplot_1.1.3